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Abstract 

The strange quark content of the nucleon (iV|ss|iV) is calculated in dynamical lattice QCD 
employing the overlap fermion formulation. For this quantity, exact chiral symmetry guaranteed 
by the Ginsparg- Wilson relation is crucial to avoid large contamination due to a possible operator 
mixing with uu + dd. Gauge configurations are generated with two dynamical flavors on a 16 3 x 32 
lattice at a lattice spacing a ~ 0.12 fm. We directly calculate the relevant three-point function 
on the lattice including a disconnected strange quark loop utilizing the techniques of the all-to-all 
quark propagator and low-mode averaging. Our result /t s =m s (N\ss\N) /M7v = 0.032(8) s t a t(22) sys , 
where m s and Mn are strange quark and nucleon masses, is in good agreement with our previous 
indirect estimate using the Feynman-Hellmann theorem. 
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I. INTRODUCTION 



In the naive quark model, the nucleon consists of three valence up and down quarks. 
This picture is made more precise by taking account of quantum effects based on quan- 
tum chromodynamics (QCD), the fundamental theory of strong interaction, with which one 
expects additional effects due to the gluon and sea quark degrees of freedom. In fact, in 
high energy hadron scatterings, these effects are observed as parton distributions of the 
gluon and sea quarks, which can be analyzed using perturbative calculations of QCD. At 
low energy, quantitative calculation of the sea quark effect is far more difficult because of 
the nonperturbative nature of QCD. In this work, we consider the nucleon strange quark 
content (iV|ss|iV). This matrix element directly measures the effect of sea quark, because 
there is no valence strange quark in the nucleon. 

The nucleon strange quark content represents the effect of strange quark on the mass of 
the nucleon, which is often parametrized by 

m s (N\ss\N) 

fT ° ~ — w N — ' (1) 

where m s and are the masses of strange quark and nucleon. ss is the scalar operator 
made of strange quark fields. This parameter is also relevant to the dark matter searches, 
as one of the candidates — neutralino in the supersymmetric models — may interact with 
the nucleon most strongly through its strange quark content via the Higgs boson exchange 
diagram IMS]. The magnitude of the matrix element (iV|ss|iV) is therefore directly related 
to the sensitivity of the present [s, 7] and future experiments. 

Another quantity of physical interest is a ratio of strange quark and light (up and down) 
quark contents: 

(N\uu + dd\N)' V ; 

The denominator (N\uu + dd\N) corresponds to the nucleon a term, which is relatively 
well-determined as it is related to an amplitude of the pion-nucleon scattering. This is 
not the case for (jV|ss|iV), for which only lattice QCD can potentially make a quantitative 
prediction. 

The numerical calculation of the matrix element (iV|ss|./V) on the lattice is not straight- 
forward, because it involves a disconnected quark-loop diagram shown in Fig. [TJ With the 
conventional method to calculate the quark propagator in lattice QCD, the computational 



FIG. 1: Disconnected three-point function relevant to (JV|ss|iV). Lines show quark propagators 
that are dressed by virtual gluons and sea quarks in QCD. The connected three lines correspond 
to the nucleon propagation and the disconnected loop arises from the strange scalar operator ss. 



cost to obtain the disconnected quark loop is prohibitively high, as one has to perform 
an expensive inversion of the Dirac operator for sources located at every lattice sites; the 
computational cost is then proportional to the lattice volume squared, {N^) 2 . Furthermore, 
since the scalar operator ss may have nonzero vacuum expectation value (VEV), which is di- 
vergent when m s is finite, one has to subtract this VEV contribution to extract the physical 
matrix element (iV|ss|./V). This requires a large cancellation that induces a large statistical 
error. 

In this work, we overcome these practical difficulties in the lattice calculation by using the 
methods of the low- mode averaging 8|, [9] and the all-to-all propagator 10, LL1[ . The all-to-all 
propagator allows us to calculate the propagation of the quark between arbitrary lattice sites 
at once, by introducing a stochastic estimator (for a practical implementation, see below). 
Although it introduces additional statistical noise, the low-mode averaging eliminates the 
noise for physically relevant low-lying quark-mode contributions and improves the statistics 
by averaging over space-time lattice sites. These techniques are crucial for the calculation 
of the disconnected diagram in lattice QCD. 

Another important advantage of this work over the previous lattice calculations of 
(iV|ss|iV) 12|dlq| is the use of a lattice fermion formulation that preserves exact chiral 
symmetry at finite lattice spacings. For both sea and valence quarks we employ the overlap 



fermion 
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171 ]. which satisfies the Ginsparg- Wilson relation [18] and thus has a symmetry 



under a modified chiral transformation 19(. This exact chiral symmetry prohibits the oper- 
ator mixing under the renormalization between ss and uu + dd, where the matrix element of 



the latter operator involves the connected diagram contribution. With the Wilson fermion 
formulation that has been used in the previous works, the operator mixing is induced due 
to the explicit chiral symmetry breaking on the lattice. Since the connected diagram con- 
tribution of uu + dd is larger than the disconnected one by an order of magnitude, this may 
give rise to a large systematic error unless the mixing contribution is subtracted nonpertur- 
batively. 

In our previous work {^(J, we used a technique to extract (N\uu + dd\N) and (N\ss\N) 
from the quark mass dependence of the nucleon mass using the Feynman-Hellman theorem. 
Since the number of sea quark mass values in the simulations was limited, the method 
had an inconsistency that the disconnected contribution was evaluated at up and down 
quark masses, which are different from the physical strange quark mass. In the present 
work, this limitation no longer remains. Although the calculation is done on two-flavor 
QCD lattices, which are available from the project of the dynamical overlap fermion by 



the JLQCD-TWQCD Collaboration 21j, an extension to the realistic 2+1-flavor QCD is 
straightforward and in fact underway. 

This paper is organized as follows. In Sec. HIl our simulation setup and the methods 
of the all-to-all propagator and the low-mode averaging are described. We investigate the 
efficiency of the low-mode averaging by comparing the statistical error of the nucleon two- 
point function as presented in Sec. IHIl Extraction of the strange quark content from the 
disconnected three-point function is discussed in Sec. [TV] Section [V] is devoted to a discussion 
of chiral extrapolation to the physical quark masses. In Sec. IVIl we emphasize an important 
role of chiral symmetry in the calculation of the strange quark content. We also make a 



comparison with previous works including the recent results 22 1 
given in Sec. IVHI A preliminary report of this work is found in 



24] 



231 ] . Our conclusions are 



II. SIMULATION DETAILS 



A. Simulation Setup 



On a four- dimensional Euclidean lattice we simulate QCD with two flavors of degenerate 
up and down quarks. As the lattice formulation, we use the Iwasaki gauge action and the 
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overlap quark action. The overlap-Dirac operator is given by 16l. Il7] 



D(m) = \ m + —j + \ m - — j 75 sgn [H w \ , (3) 

where Hy/ = 7 5 D^(— mo) is the Hermitian Wilson-Dirac operator and mo = 1.6 in this 
study. The mass parameter m corresponds to the up-down or strange quark mass. We 
also introduce an additional Boltzmann factor 25| which does not change the continuum 



limit of the theory but substantially reduces the computational cost to calculate sgn[if w ] by 
prohibiting the exact zero modes and suppressing near-zero modes of Hy/. This additional 
Boltzmann factor induces a side effect that the global topological charge Q during the hybrid 
Monte Carlo update is fixed. We simulate only the trivial topological sector Q = in this 
study; the effect of fixing topology is suppressed by an inverse power of the space-time 



volume l/(NgN t ) 26j and turns out to be small (typically below a few percent level) in our 



studies of meson observables 
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29l | . We expect that it is even smaller for baryons. 
Our gauge configurations are generated on a N% x N t = 16 3 x 32 lattice at a gauge 
coupling (3 = 2.30 where the lattice spacing is determined as a = 0.118(2) fm using the 
Sommer scale ro = 0.49 fm as an input. We accumulate 100 independent configurations of 
two-flavor QCD at three values of up and down quark masses m u d = 0.025, 0.035, and 0.050, 
which cover a range of the pion mass M n = 370-520 MeV. The physical quark masses are 
fixed as m u d tP hys = 0.0034 and m x r hy S = 0.077 from our analysis of the pion and kaon masses 
We refer the readers to [21f| for further details of the configuration generation. 
We take two values of the valence strange quark mass m S)Va i = 0.070 and 0.100 close to 
^s, P hys, and calculate two- and three-point functions 

C 2 r pt (y,t src ,At) = j!pJ2 tTs [r(iV(x,t src + At)iV(y,t src ))] , (4) 

S X 

C 3 r pt (y, t src , At, At,) = Yl ( tr * [ r ( iV ( x ' *«c + At)S lat (z, t src + At,)N(y, t s 



x,z 



-(S lat (z,t src + Ai s ))tr s [r(iV(x,t src + At)iV(y,t sl . c ))]}, (5) 



where we use the nucleon interpolating field N = e abc {u^C^db)u c with the charge conjuga- 
tion matrix C = 7472- The trace "tr s " is over spinor index of the valence nucleon and (■ • • ) 
represents a Monte Carlo average. The scalar operator made of the strange quark field is 
given by 



on the lattice for the overlap-Dirac operator ([3]). To obtain the continuum operator S cont (iJ,) 
at the energy scale /i, we need the renormalization factor Zs(p) as S cont (fi) = Zs(fi)S la,t . 
The details including possible operator mixing are discussed in Sec. ED 

We take two choices of the projection operator F — F± — (1 ± 74)/2, which correspond 
to the forward and backward propagating nucleons, respectively. The two- and three-point 
functions are averaged over the two choices of T 

C 2pt (y, t 8IC , At) = - {c 2 r + (y, t SIC , At) + C 2 r pt (y, t src , N t - At) } (7) 
C 3pt (y,t src ,At,At s ) = ^{cl+{y,t SIC ,At,At s ) + Cl pt {y,t src ,N t -At,N t -At s )} (8) 

in order to reduce statistical errors. 



B. All-to-all quark propagator 

The three-point correlation function C 3pt is calculated by appropriately connecting the 
quark propagator as shown in Fig. [TJ The conventional method to calculate the 

quark propagator is not suitable to construct the disconnected quark loop starting from and 
ending at arbitrary lattice sites since the source point y has to be fixed at a certain lattice 
site. Indeed, we use the all-to-all quark propagator technique, which enables propagations 
from any lattice site to any site, following the strategy proposed in [10, 1 1 it ] - 

It is expected that low- lying eigenmodes of D(m) dominantly contribute to the low- 
energy dynamics of QCD. We calculate the low-lying eigenvalues and eigenvectors using the 
implicitly restarted Lanczos algorithm, from which we can construct their contribution to 
the quark propagator exactly as 

N e 1 

(D-\m)) low (x,y) = £ j^v^i^ivV, (9) 
i=i ^ ' 

where A^(m) and v®(x) represent the i-th lowest eigenvalue and its associated eigenvector 
of D(m), respectively. Note that the eigenvectors are independent of valence quark masses. 
The number of low-lying eigenmodes N e we calculated is 100 in this study. 

The remaining high-mode contribution is estimated stochastically. We prepare a single 
Z 2 noise vector r](x) for each configuration and split it into = 3 x 4 x N t /2 vectors r]^ d \x) 
(d = l,...,Nd), which have nonzero elements only for a single combination of color and 
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spinor indices on two consecutive time slices. The high-mode contribution is then estimated 

as 

N d 

(D-\m)) high (x,y) = £l>V>(x)r,to(yy, (10) 

d=l 

where ift( d '(x) is obtained by solving a linear equation for each noise vector 

D{m)^ (d \x) = (1 - V low )r] id \x) (d = 1, N d ). (11) 
Pi ow is a projector to the subspace spanned by the low-modes 

V low (x,y) = J2v {r) (x)vV( y y. (12) 

i=l 

We use this all-to-all propagator, namely, (jUJ) plus ffTOj) . to calculate the disconnected quark 
loop and the vacuum expectation value of S' lat in C^. 

C. Low- mode averaging 

In principle, we can use the all-to-all propagator to calculate nucleon correlators, namely, 
C 2p t and the piece representing the nucleon propagation in C 3pt . However, these quantities 
decay exponentially as the temporal separation At increases, so that the contributions to 
the nucleon correlator from the high-modes (1101) are not sufficiently precise at large At when 
we take only one noise sample for each configuration. 

In this study, we therefore use the low-mode averaging (LMA) technique proposed in 
jsl, Suppose that we decompose the conventional quark propagator into its low-mode 
part, which is in the subspace spanned by the low- modes and the remaining high- mode 
part. We can then write C^pt in terms of the following eight contributions: 

<-"2pt — ^2pt ' 2pt ' 2pt "+* 2pt "+* 2pt + 2pt + 2pt "+" °2pt ■ 1 16 ) 

Here, is constructed only by the low-mode part of the quark propagator; C l ^ t is the 
one in which two of the valence quarks are made of low-modes and the other is the high- 
mode part. The other combinations are understood in a similar manner. Since the ensemble 
average can be taken for each term of (fT51) . we attempt to reduce the statistical error for 
individual contributions. 

Relying on the translational invariance, we may replace by a more precise estimate 
by averaging over the location of the nucleon source point (y, £ src ). No additional inversion of 



the Dirac operator is necessary to take the average, as we can explicitly use the representation 
(Q made of low-mode eigenvectors. This LMA technique is very effective in reducing the 
statistical error of C 2pt at large At when C 2pt is well dominated by C^. 

In this study, we employ LMA to calculate C 2pt and the nucleon piece of C^ p t- We also 
test an extension in which additional three contributions, C^, C^, and C^, are averaged 
over the source location by using the all-to-all propagator. The signal may be improved if 
the reduction of the statistical error by the source average outweighs the induced noise from 
the high-modes. The result of the test is shown in the next section. 

D. Smeared nucleon operators 

Since C 2pt and C^pt decay quickly as a function of At, we need to use smeared nucleon 
operator that suppresses excited-state contaminations at small At. 

For the (local or smeared) quark field, we consider the following three choices: 

1. local 

ffloc(x,t) = ?(x,t). (14) 

2. exponential smearing 

C p r (x,t) = ^exp(-B|r|) g (x + r,t), (15) 

r 

where the parameter B is set to 0.350, 0.375, 0.400 at m ud = 0.025, 0.035, 0.050, 
respectively. 

3. Gaussian smearing 

= W (l + ^hY) q(y,t), # x , y = E(<W-; + <Wh), ( 16 ) 

y ^ > x,y j = i 

where the parameters u> = 20 and iV = 400 are chosen so that the extent of the 
smeared operator is roughly equal to that of ( fl5i) with B = 0.400. 

Then, the nucleon interpolating fields A^ oc (x, t), i\r^(x, t), A^ s r (x, t), are constructed from 
the corresponding local or smeared quark fields. 

When we smear the quark field, we fix the gauge to the Coulomb gauge. With this choice 
one can avoid significant statistical noise coming from the fluctuation of the gauge link. 
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At 

FIG. 2: Effective mass Mjv(At) from the nucleoli two-point function C2 P t at m u d = 0.050. The 
local operator is used for both source and sink. Circles show the result of the conventional point 
source, while squares (triangles) are obtained by averaging the Cfpt (C^pt + ^2pt + Qjpt + ^2 P t) 
contributions. Circles and triangles are slightly shifted in the horizontal direction for clarity. 

The Gaussian smearing is particularly useful for the sink smearing, since the number of 
numerical operation ~ N x N% is smaller than ~ N® for the case of ( TT5l) . 



III. IMPROVING THE NUCLEON TWO-POINT FUNCTION 

Since the disconnected three-point function C^t is extremely noisy, it is crucial to reduce 
the statistical noise and to extract the signal at relatively small time separations. We 
therefore tested various methods to improve the signal on the nucleon two-point functions 
C 2 pt before applying them to the three-point functions. 



A. Low-mode averaging 

As mentioned in the previous section, we consider two options: (i) to average only C^L 
over the source locations, (ii) to average also C^h + C l J]i + CtfE. The second choice requires 
the high-mode of the quark propagator (-D _1 )hi g h(a ; , y), which is calculated stochastically as 
in (I10j). 

In Fig. [2] we plot the nucleon effective mass Mjy(At) with the local source and sink 
operators at our heaviest quark mass m u d = 0.050. The data without LMA (circles) show 
a rapidly growing statistical error as At increases, so that the error at At = 10 where 
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FIG. 3: Effective mass M/v(Ai) from C2 P t with an exponentially smeared source at m u d = 0.050. 
The symbols are the same as in Fig. [2j 

the plateau is approximately reached is already as large as 4%. By averaging over the 
source locations for Cfit (squares), the statistical error is reduced by a factor of about 3. 
Further improvement of a factor of 2 is possible if we average over the source points also for 
C"pt + ^2pt + ^2pt' as shown by triangles. 

A similar comparison of Mat (At) at m u d = 0.050 but with the exponentially smeared 
source and a local sink is shown in Fig. [3j (But LMA is done over a limited number of the 



source location iV sl 



N f x 16. For discussions, see below.) We observe that LMA for C. 



in 

»src — -» v t •*> J-U. J-'ui Luauusaiuiis, sec ueiuw.j we uuaei ve inai uivii-v iui ^2pt 

is efficient when combined with the smeared source, while the effect of the extended LMA 
for C*2p t + Cgpt + CM is not substantial, i.e., the reduction of statistical error is only about 
30%. 

Although the effect of LMA to reduce the statistical noise is significant, it is also true 
that it requires substantial computational effort. If we average over the entire space-time 
source points, the computational cost scales as (N^ x N t ) 2 , which is prohibitive unless we 
use the fast Fourier transform. If we combine LMA with the smeared source, another factor 
of is necessary, which is not feasible any more. We therefore consider averaging over a 
limited number of source locations. Since the correlators from different source points are 
statistically highly correlated, this might not spoil the efficiency of LMA largely. 

In Fig. HI we compare the data of Mjy(At) obtained using LMA with a different number 
of source points averaged iV SI . c . The plot shows the results of LMA for both C|p t and 
^22pt+^2pt+^2pt w ith iV src = N t (squares), N t x8 (triangles down), and N t xl6 (triangles up). 
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FIG. 4: Comparison of M^(At) obtained with different numbers of source locations for LMA. 
Circles are those without LMA. Results averaged over the time slices are shown by squares. We 
obtain down- and up-triangles by further averaging over 8 and 16 spatial sites at each time slice, 
respectively. In the plot, iV src represents the number of the source locations. 

For N src = N t , the spatial location of the source is fixed and the average is taken over N t time 
slices. For N STC = N t x 8, points of spatial coordinates or N s /2 in three spatial dimensions 
are all averaged; for N t x 16, we also average over (iV s /4, N s /A, N s /A), (N s /4:, N s /A,3N s /i), 
(Ay4,3iV s /4,3Ay4), and (3iV s /4, 3iV s /4, 3iV s /4) (and all possible permutations) for each 
time slice. 

From Fig. 4 we observe that the result with N STC = N t is already very good, while the 
improvement with N SIC — N t x 8 is not substantial. Beyond this number, we do not gain 
significant improvement. Note that the maximal number of points we took iV src = Nt x 16 
corresponds to the data shown in Fig. [3] (triangles). 

Overall, taking the cost of numerical calculation into account, the best choice would be 
N STC ~ N t x 8; in our following analysis we choose A^ src ~ N t x 16, which has been still 
doable. The advantage of LMA for C^L is always clear, while that for + C l 2 p[ + C%£ t 
depends on the channel or source operator. Therefore, we average only C%L when we use 
the smeared sink, which is numerically more costly. 
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FIG. 5: Effective mass Mjv(At) with the Gaussian smeared source and sink at m u d = 0.050. Circle 
are without LMA; triangles are obtained by averaging over 16 spatial sites at each time slice. 

B. Sink smearing 

The smearing of the source operator is routinely used in many lattice calculations. It is 
designed to deplete the overlap with excited-state contributions so that the plateau of the 
effective mass constructed from the two-point correlator appears earlier in At. By using 
the smeared operator also for the sink we expect that the excited-state contaminations are 
further reduced, but usually the benefit is not clearly seen mainly because the statistical 
noise increases with the smeared sink. Since the numerical cost for the sink smearing is high 
in general [~ (iVf) 2 ], it has not been commonly used. 

The situation may be different for three-point functions, where an operator is inserted in 
the middle of the two-point function. Here the nucleon and its excited states are created at 
the smeared source point and propagate until the point of the operator is reached. Between 
these two points, the depletion of the excited states is at work because of the smeared source. 
After the insertion of the operator, the nucleon and its excited states propagate until they 
are absorbed by the sink. In this second propagation, the excited states are not necessarily 
suppressed, since the operator insertion may excite the nucleon, i.e., (N\ss\N') ^ 0, and the 
sink operator may have substantial overlap with the excited state \N'). This is indeed the 
case in our calculation of the three-point function relevant to the strange quark content, as 
we will see in the next section. 

We therefore utilize the smeared operator also for the sink. Since the conventional choice 
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FIG. 6: Improving the statistics by averaging the nucleon two-point functions calculated from 
four different source points at time slices t src = 0, 8, 16, and 24. The result for Mjv(Ai) (squares) 
is compared with that without the average, i.e., t src = 0. The quark mass is m u d = 0.050. The 
Gaussian smearing is used for both the nucleon source and sink. 

?srar( x >^) CGD requires a numerical cost proportional to N% for each (x, t), we use gf^ r (x, t) 
(TIE]) , instead. Figure [5] shows Mjy(A£) with this Gaussian smeared operator for both the 
source and sink. Although the statistical signal is worse compared to the case of the local 
sink shown in Fig. [3] and HI we may improve it using LMA for as shown in Fig. [5] by 
triangles. Further improvement is not expected with the average over Cfj^ + C l 2 ^ t + C^, as 
in the case of the smeared source and local sink (Fig. [3]). 



C. Duplication 

Instead, we simply repeat the calculation 4 times by setting the source at different time 
slices. Namely, we calculate the nucleon two-point function locating the source on the time 
slices t src = 8, 16, and 24, in addition to the original choice t src = 0, and average over 
these duplicated correlators. The effect is shown in Fig. El where we observe a reduction 
of the statistical error by a factor of 2 at large time separations. However, we find that 
the further average of the duplicated correlators is not substantial. This is tested at m ud = 
0.025 by calculating the nucleon two-point function locating the source on the time slices 
t src = 4, 12, 20, and 28 besides t BIC = 0, 8, 16, 24. Therefore, we restrict the number of the 
duplication of the nucleon two-point function for other quark masses. 
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FIG. 7: Noise-to-signal ratio of the nucleon correlator at m u d = 0.015, 0.025, 0.035, and 0.050. 
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TABLE I: Choices of the scheme of averaging the nucleon correlator in this work. For different 
smearing operators at the source and sink, we list the number of source points iV src averaged in 
LMA, the contributions to the correlator averaged in LMA (C^ t or C^ t + C^h + C%g t + C^ t ), 
and the number of the duplications of the conventional correlators. 

Figure [7] shows the increase of the statistical noise in C*2 P t for the case of smeared source 
and sink. The plot shows the data at four different quark masses m u d = 0.015, 0.025, 0.035, 
and 0.050. As expected, the noise grows more rapidly for lighter quarks. Since the plateau in 
the effective mass is reached at around At = 5, we need at least At = 10 in the calculation 
of the three-point functions. At the lightest quark mass m u d = 0.015, the error around 
At = 10 is too large (~ 10%) to be useful in the analysis of the disconnected three-point 
functions. We therefore discard this data point in the analysis of the strange quark content. 

In order to optimize the statistical signal in the calculation of the disconnected three- 
point function for a given amount of computer time, we choose different schemes of averaging 
the correlators depending on the source and sink smearing combinations. These include the 
choices of the contributions averaged in LMA (C%L or C l 2 l L + C^A + C& t + C^pt), the number 
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FIG. 8: Ratio R(At,At s ) with At = 11 at m ud = 0.050 and m SjVa i = 0.100. Circles (triangles) 
are results obtained with the local (exponentially smeared) source and local sink, whereas squares 
are calculated using the Gaussian smeared source and sink. The vertical lines show the locations 
of the nucleon operators. The noisy high-mode contribution to the quark loop is ignored in this 
plot. 

of source points N STC averaged in LMA, as well as the number of the duplications of the 
conventional correlators. Our choices in this work are listed in Table [B 



IV. EXTRACTION OF THE STRANGE QUARK CONTENT 
A. Finding a plateau in the three-point function 

We extract the strange quark content on the lattice (A r |S' lat |A r ) from a ratio of 
C 3pt (At,At s ) andC 2pt (At) 

R(At,At s ) = ^(At At,) y {N \s^\ N ) (17) 

C2 p t(^t) At,At s ->oo 

where At is the temporal interval between the nucleon source and sink. The scalar operator 
5* lat is set on the time slice apart from the nucleon source by At s . Note that C3 P t(At, At s ) 
and C 2pt (At) are calculated with LMA. We suppress the coordinates of the nucleon source 
location (y,t S rc) presented in ([7]) and (jS]). 

In order to extract (7V| *S' lat | ^V) , we first have to identify a plateau in the ratio R(At, At s ) at 
sufficiently large At and At s . For this purpose, we look at the same ratio but approximated 
by taking only the low-mode contribution in the strange quark loop. Namely, the piece 
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FIG. 9: Same as Fig. [8j but for m u d = 0.035 (top panel) and m uc i = 0.025 (bottom panel). 



of S lat (z) in (jSJ) is replaced by its low-mode contribution Tr[(D _1 (m))i ow (2;, z)\. We expect 
that the ratio R(At, At s ) is dominated by this low-mode contribution, because the high- 
mode contribution that leads to the ultraviolet divergence in the continuum limit cancels 
by the VEV subtraction in ([5]). Low-energy physics must be well described by the low- 
mode contribution in the strange quark loop. This approximation is finally removed in our 
calculation by the full calculation, but here we consider the approximately calculated ratio 
to identify the plateau, where the ground-state nucleon dominates. 

Figure [8] shows the approximated ratio obtained at m ul i = 0.050 and m SyVa i = 0.100 with 
various combinations of the source and sink smearing. The separation between the source 
and sink is fixed to At = 11, and the location of the scalar operator At s is varied. Thus, we 
expect a signal around At s ~ At/2. We observe a plateau between At s = 3 and 8, when the 
source and sink operators are both smeared with the Gaussian smearing (jlBI) . as shown by 
filled squares. The data with the local source and sink (open circles) show a slight increase 
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• w/ high-mode 
□ w/o high-mode 



At 

FIG. 10: Ratio R(At = 13,Ai a ) at m ud = 0.050 and m S)Va i = 0.100 with (filled circles) and 
without the high-mode contribution to the strange quark loop (open squares). The horizontal lines 
show the result of a constant fit R(At) and its error band. 

in the same region but do not reach the value of the plateau for the smeared source-sink 
combination. 

The data of the smeared source and local sink (open triangles) show a bump around 
At s ~ 2 — 6 and decrease towards At s = 11, so that the plot looks asymmetric. This can 
be explained by an excited-state contamination on the sink side (At s = 11) because the 
sink operator is local. Therefore, unlike the case for the two-point function, the use of the 
smeared operator for both source and sink is essential for the three-point function in order 
to extract the ground-state signal. 

Similar plots are shown for m U d = 0.035 and 0.025 in Fig. [9j We observe similar behavior 
of the approximated ratio. 



B. Bare results for the strange quark content 

The ratio R(At, At s ) in (fT7|) without the low-mode approximation is shown in Fig. [10] 
(filled circles) together with that of the low-mode approximation (open squares). Here, 
the data for At = 13 are shown. Although the statistical noise is much larger when the 
high-mode contributions are included, the central value is unchanged. 

Since the high-mode contributions are calculated with random noise (Tl0|) . the larger noise 
is expected. But, because the noise given for each time slice is statistically independent, the 
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m =0.050,m =0.100 



10 



12 



FIG. 11: Results of the constant fit for R(At, At s ) in the range At s = [5, At — 5]. The data at 
m ud = 0.050. The left and right panels show those at m s ^ va i = 0.070 and 0.100, respectively. 

correlation among the data points at different At s is expected to come mainly from the low 
modes, provided that the high-mode contribution to the ratio is negligible, which is indeed 
the case within our statistical accuracy. The statistical error is then effectively reduced by 
averaging over different At s . In Fig. [IUJ the result of a constant fit for At s = [5,8] is shown 
by a horizontal line together with a band showing the resulting statistical error. In this case, 
the statistical error of the fitted value is about a half of that of each point, because four 
data points are averaged. We also checked that the statistical correlation among the points 
at different At s is an order of magnitude smaller than the variance of each point. 

For the final result, we take the full data including the high modes and fit in the region 
where the approximated ratio shows a plateau. To be specific, we fit in the region At s = 
[5, At -5] with At > 11. 

Figures [TTHT3] show the results of the constant fit for each At. We find that the results 
are stable under the change of At. We then fit these results by a constant in At = [11,15]. 
The statistical error is estimated using the jackknife method. The numerical results are 
listed in Table HH 

In order to estimate the systematic effect due to possible contamination of the excited 
states, we also test a fitting form for R(At, At s ) taking account of the first excited state: 



R(At, At s 



c - Cie 



-{2M +AM)At/2 



cosh(AM(At s - At/2)), 



(18) 



where the first and second terms represent the contributions from the ground and first 
excited states, respectively. AM is the mass gap between these two states. To make this 
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FIG. 12: Same as Fig. [IT] but at m ud = 0.035. 
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FIG. 13: Same as Fig. OH but at m ud = 0.025. 



fit stable, we carry out a simultaneous fit in terms of At s and At using a slightly wider fit 
range, At s = [4, At — 4] and At > 11. We also use the ground-state mass Mo determined 
from the nucleon two-point function. The excited-state contribution represented by the c\ 
term turned out to be small: in the maximum case (At = 11) it is about 0.04(8) compared 
to the main contribution cq ~ 0.3(1). For large At, the excited-state contribution is more 
suppressed. This is expected from the small At s and At dependence of the ratio shown in 
Figs JTOHTBl We therefore use the results in Table HT1 in the following analysis without adding 
further errors due to the excited states. 
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mud Fit range of At m S;Va i = 0.070 m S)Va i = 0.100 
0.050 [11,15] 0.345(89) 0.286(83) 
0.035 [11,15] 0.089(100) 0.070(96) 
0.025 [11,15] 0.351(128) 0.303(126) 



TABLE II: Strange quark content (A r |S' lat |A^) calculated on the lattice at each quark mass. The 
fit range of At is also listed. 
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FIG. 14: The dependence of (A r |5 lat |A^) on the up and down quark mass m u d (given in the lattice 
unit). Open circles and triangles are the data at each m u d and m s = 0.070 (circles) and 0.100 
(triangles). The data linearly interpolated to the physical strange quark mass m s . p h ys is shown by 
filled diamonds. Dashed and solid lines show the fit curve at m SiP hy S obtained from the constant 
and linear extrapolations. 

V. CHIRAL EXTRAPOLATION TO THE PHYSICAL POINT 

In this section, we discuss on the extrapolation of our lattice data to the physical quark 
masses. We have three data points corresponding to up and down quark masses m u d in 
the range of M w = 370-520 MeV. For the strange quark mass we have two data points 
sandwiching the physical strange quark mass. 

Our data for the matrix element (N\S \N) are plotted as a function of in Fig. [TH 
We do not observe statistically significant dependence of (iV|S' lat |iV) on both m u d and m s . 
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X 2 /d.oi. 


d.o.f. 


co 




Cl,s 


(A|<S lat |iV) 


constant 


1.63 


5 


0.24(6) 






0.24(6) 


linear 


2.39 


3 


0.22(24) 


3.5(5.7) 


-1.44(52) 


0.12(22) 



TABLE III: Numerical results of chiral extrapolation. We also list (iV| *S' lat | N) extrapolated to the 
physical point. 

By fitting the data linearly in m u d and m s as 

(N\S lat \N) = c + c 1M m ud + c ltS m s , val , (19) 

we obtain the numerical results of the fit parameters Cq, ci )1W j, and C\ yS listed in Table IHIl 
We also show the result of a constant fit including only the Cq term in (fl~9]) . Both results are 
consistent with each other, but the linear extrapolation gives a larger error at the physical 
point. 

Assuming that the quark mass dependence of the nucleon mass is reliably described by the 
chiral perturbation theory, we also attempt an extrapolation using the formula provided by 
the SU(3) heavy baryon chiral perturbation theory (HBChPT). From the chiral expansion 
of Mjv 3jj and the Feynman-Hellmann theorem (|3T|) . which will be discussed in Sec. IVIl 
the quark mass dependence of (A r |S' lat |A r ) up to the next-to- leading order is given by 

(iV|S lat |iV) = -c„-B !^C NNK M K + 2C NNT] M T] J , (20) 

where the coefficients Cnnk and Cnnt) are written as 

1 (5D 2 -6DF + 9F 2 ) 
C NNK - — 2 , (21) 

r 1 — ~ 3F)2 (22) 

C NNv - — 2 . (22) 

The axial couplings F and D are phenomenologically well determined and we fix them as 
D = 0.81 and F = 0.47 32]. For the pseudoscalar meson masses Mk and M v , we use 
the Gell-Mann, Oakes, and Renner (GMOR) relations M\ = B(m U( i + m s ) and M 2 = 
2B(m u d + 2m s )/3, which are valid at the leading order of the quark masses. We fix the 
low-energy constants / and B to the values obtained in our study of the pion mass and 



decay constant 28j. Note that the contributions of the decuplet baryons are ignored in this 
analysis. 
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FIG. 15: The chiral fit of (JV|5 lat |./V) based on the next-to-leading order HBChPT ((20]). Solid and 
dashed lines show the fits at m StVa i = 0.070 and 0.100. 



1 ' 


> 


1 1 

C 


) 


i '. 


1 o m sv a] =0 - 070 

: A m svaI =0.100 

I 









■ 4 I i I i I i I i I i I i I 

' 0.01 0.02 0.03 0.04 0.05 0.06 



FIG. 16: The chiral fit using f|23|) with a higher order term. 

As one can see from Fig. [15], this function does not describe the numerical data; the value 
of x 2 P er degree of freedom (d.o.f.) is unacceptable (~ 20). The main reason is that there 
is no free parameter to control the quark mass dependence, i.e., the coefficients of Mk and 
M v in ( 120]) are completely determined phenomenologically. In other words, if we leave / as 
a free parameter for instance, the resulting value is unreasonably large. 

If we add a higher order analytic term as 

(iV|£ lat |iV) = -c s -B !^C NNK M K + 2C NNn M^ + c 2 M 2 K , (23) 

the fit becomes reasonable as shown in Fig. [16] for which x 2 /d.o.f. is acceptable (~ 1.9). 
Fit parameters obtained with ( 120]) and ([23]) are summarized in Table IIVI The resulting fit 
parameters suggest that the chiral expansion does not converge well. In fact, if we look at 
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X 2 /d.o.f. 


d.o.f 


-C s 


C2 


(N\S lat \N) 


Equation (|20j) 


19.5 


5 


5.48(6) 




1.24(6) 


Equation ([23]) 


1.88 


4 


2.82(23) 


21.2(1.8) 


0.28(10) 



TABLE IV: Numerical results of chiral fits using the SU(3) HBChPT formulas, i.e., (HOD and ([23]). 

the individual contributions to (A r |S' lat |A^) from each term in ( [23]) . all of them are an order 
of magnitude larger than the data themselves, and the final result is obtained by a large 
cancellation. 

Because of this poor convergence of the chiral expansion, we use the result of the HBChPT 
analysis only to estimate the systematic uncertainty. Namely, we take the result from the 
constant fit as a central value of (A r |S' lat |A r ) at the physical quark masses. The systematic 
error due to the chiral extrapolation is estimated by a difference from the results of the linear 
([19]) and HBChPT fits ([23]). Then, we obtain (iV|,S lat |iV) = 0.24(6)(16) at the physical quark 
masses. The first and second errors represent the statistical and systematic ones. 

Using the experimental value of M^, this is converted to the strange quark mass contri- 
bution to M N defined in ([1]) as 

f Ts = 0.032(8)(22). (24) 

Since the combination m s S' lat is invariant under renormalization, no renormalization factor 
is required to obtain ( [241) . 

The y parameter ([2]) is defined as a ratio of the strange and ud quark contents. We obtain 

y = 0.050(12) (34), (25) 

where we use an estimate (N\uu + dd\N) = 9.40(41) for the denominator, which is taken 
from our study of the nucleon sigma term [20|. 

A simple order counting suggests that the discretization effect is 0((aA) 2 ) ~ 9% when we 
take A ~ 500 MeV. Other systematic errors including those of finite volume effects would 
not be significant, given that the statistical and systematic errors in ( ]24l) and ( [25]) are so 
large (~ 70%). 
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FIG. 17: The disconnected diagram contributing to the renormalization of flavor-singlet scalar 
operator (cross). At higher orders, the quark loop and the quark line on the bottom may be 
connected by an arbitrary number of gluon propagators. Since the quark-quark-gluon vertex con- 
serves chirality, the chirality of the quark propagating in the loop does not change, as far as the 
regularization respects chiral symmetry. 

VI. COMPARISON WITH PREVIOUS LATTICE CALCULATIONS 

In this section, we emphasize an important role played by the exact chiral symmetry in 
the calculation of the strange quark content. Then we compare our result with the previous 
calculations. 

A. Renormalization issue of the operator ss 

First, let us consider the renormalization of the ss operator in the flavor 577(3) symmetric 
limit for simplicity. Using the flavor triplet quark field ip, the ss operator can be written in 
terms of flavor- singlet and octet operators as 



where A 8 is a Gell-Mann matrix. Note that, in this section, we put the superscript "phys" 
on the renormalized quantities denned in the continuum theory to distinguish them from 
bare operators, which is in our case defined on the lattice. 

In general, the singlet and octet operators may be renormalized differently 




(26) 



(#) phys = Zo{W>), 



(27) 
(28) 
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with different renormalization factors Z$ and Z%. Here, we assume that the chiral symmetry 
is preserved in the renormalization scheme used to calculate (|27|) and f l28|) . Otherwise, there 
is a mixing with lower dimensional operators for the flavor-singlet operator ( 12 7p , as discussed 
below. The operator (ss) phys is then expressed in terms of bare operators as 

(ss) phys = \ {(Z + 2Z 8 )(ss) + (Z Q - Z 8 )(uu + dd)} , (29) 

which implies that the ss can mix with uu + dd unless Z$ = Z 8 . The difference Z — Z 8 
arises from disconnected diagrams such as those shown in Fig. [TTJ, which exist only for the 
flavor- singlet operator. 

When the renormalization scheme respects chiral symmetry, the disconnected diagrams 
vanish in the massless limit, because the quark loop starting from and ending at a scalar 
operator ss = slSr+srSl has to change the chirality in the loop while the change of chirality 
does not occur by attaching any number of gluon lines to the quark loop. It means that 
Z = Z 8 is satisfied for mass independent renormalization schemes, as far as they maintain 
exact chiral symmetry. This also applies in the case of the overlap fermion formulation on 
ihe lattice, as there is an exact chiral symmetry guaranteed by the Ginsparg- Wilson relation 
18] at finite lattice spacings jl^l ]. 

Thus, the renormalization of the scalar operator reduces to a multiplicative renormaliza- 
tion (ss) phys (n) = Zs(n)S lat with Z$ = Z = Z 8 . Here we specify the renormalization point 
/x for the renormalized operator (ss) phys . The value of Zs{fi) is nonperturbatively calculated 
in {33} as Zs(2 GeV) = 1.243(15) on our lattice. For the numerical results of fx a (124")) and y 
( )25l) quoted in the previous section, the renormalization factor is unnecessary, because they 
are related to a renormalization invariant operator m s ss or a ratio ss/(uu + dd). 

As it is clear from the above discussion, the explicit violation of chiral symmetry with 
the conventional Wilson-type fermions induces a mixing between the strange and ud quark 
contents. In addition, the flavor- singlet scalar operator mixes with an identity operator, so 
that (1291) is modified as 



[ss) phys - - 



[Z + 2Z 8 )(ss) + (Z Q - Z 8 )(uu + dd) + h + 

a 6 



(30) 



3 

where the term bo /a 3 represents the power divergent mixing contribution. This contribution 
from the identity operator must be subtracted as a part of the vacuum expectation value of 
ss. Because of the cubic divergence, this results in a large cancellation toward the continuum 
limit. 
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Furthermore, since Z$ — Z% does not vanish when chiral symmetry is violated, ss mixes 
with uu + dd, which induces a connected diagram contribution in the calculation of the three- 
point function. Since the connected diagram is larger than the disconnected contribution by 
an order of magnitude, the whole effect from (Z — Z s )(uu + dd) could be substantial, even 
though the difference Z — Z 8 may be small. This possibility has been neglected in most of 
the previous lattice calculations using the Wilson-type fermions. 



B. Direct and indirect calculations 



The strange quark content can also be calculated from the m s dependence of M N through 
the Feynman-Hellmann theorem 

(NMN) = ^. (31) 

We refer to this method as the spectrum method in the following. Exact chiral symmetry 
plays a crucial role in this method, too. With the explicit chiral symmetry violation, masses 
of sea and valence quarks, m/ jSea and m/ jVa i (where / distinguishes the quark flavors ud 
and s), depend on the sea strange quark mass m S)Sea . Namely, there is an additive mass 
renormalization Am 

m £c S a = 2m( m /,sca + Am), (32) 

m fZ = Z m (m /jVal + Am), (33) 

when we relate the bare quark masses on the lattice (m/ iSea and m/ iVa i) to their counter- 
parts (my 1 ^ and m^ 1 ^) defined in some continuum renormalization scheme. Z m is the 
multiplicative renormalization factor. With dynamical Wilson fermions, this additive mass 
renormalization Am is of the cutoff order, ~ 1/a, and its dependence on the sea quark mass 
is a quantity of order unity. 

Then, we can write the relevant partial derivative dM^ /<9m S)Sca calculated on the lattice 
in terms of the "physical" quark mass dependence of Mn as 

8M N _ dmgSSL dM N | d<£ s ea dM N | dmfj^ dM N 



dm StSea dm s , sea «9m? h s S <9m s , sea dm^ dm s , sea dmf/^ 



— Z m 



(iV|s S |iV) phys + 7^-——— (N | uu + dd + ss\N} phys 



(34) 
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where the matrix elements appearing on the right-hand side are those with the continuum 
renormalization scheme. The last term must be subtracted from (l/Z m )<9Mjv /dm SySca to 
obtain the strange quark content. It requires a calculation of the light quark content {N\uu + 
dd\N), which is dominated by the connected diagram, and of the m SjSea dependence of Am, 
which strongly depends on the details of the lattice action used in the calculation. In the 
literature, this subtraction was considered only in [15j, where the subtraction induced a 
rather large statistical error. 

One may avoid this problem by differentiating in terms of pion and kaon mass 
squared, M% and M^, instead of m s , assuming the GMOR relations Ml = 2Bm ud: M 2 K = 
B(m U( i+m s ). Since the quark masses appearing in the right-hand side of the GMOR relations 
contain the additive mass renormalization Am, the above subtraction is not necessary. But 
the method introduces another uncertainty, because the GMOR relations are valid only at 
the leading order of the quark mass, and the higher order terms are not negligible in general. 
This method has been applied in the analysis of [22]. 

In Fig. [18] (top panel) we compare our result ( 125]) for the y paramet er plot ted by a 



solid circle with those from previous studies u sing the Wilson-type actions 12l-ll5|. Among 



these, 



12 



131 ] are quenched calculations and 14j, ll5j contain the effects of two dynamical 
flavors. Rather large valu es y = 0.4-0.8 were obtained in the calculations from the nucleon 



three-point functions 



12 



2Hl4j, for which the above mentioned contamination was not taken 
into account and large systematic error is expected. An exception is the UKQCD's calcu- 
lation with the spectrum method [15]; the subtraction of the contamination led to a large 
uncertainty in y. 

In the same figure, we also compare our result ( 125]) with our previous estimate y = 
0.030(16) (ig) from the spectrum method 20|, where the first and second errors are statistical 
and systematic, respectively. Because of the exact chiral symmetry satisfied in both of our 
calculations, these two points are free from the contamination and consistent with each 
other. 



Recently there have been two calculations published 22|, |23j. The analysis of Young 



and Thomas 



22j fits the data from recent calculations of the baryon spectrum done by 



the LHPC [34j and PACS-CS [35| Collaborations, and takes a derivative in terms of M% 
and M\. As already mentioned, the problem of the operator mixing is avoided in this 
method and the authors obtained a result consistent with ours. Toussaint and Freeman 
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This work 



Young & Thomas (09) 
Toussaint & Freeman (09) 



1 



[A 



I A 1 



I B 1 



B 



a 



0.5 
y-parameter 



FIG. 18: (Top panel) The comparison of the y parameter with 



this study (|25p and our estimate using the spectrum method 



previous studies. The result in 



201 ] are plotted by filled and open 



12 



circles. Triangles represent previous direct calculations from the nucleon three-point functions 

Squares are the results of that used the spectrum method. We plot two results obtained 
with and without subtracting the contamination due to the chiral symmetry breaking. (Bottom 



otted. Values are converted to the y parameter using 



panel) Two recent results [22j, are p 
the quark mass ratio m s /m u d = 27.4(4) 



371 ] and the nucleon a term I20II . The quoted errors are 



statistical only except for our studies and the two recent calculations 




uses the data for the nucleon mass obtained by the MILC Collaboration using the so-called 
"asqtad" quark action, which is a variant of the staggered fermion formulation. They use a 
clever idea of extracting the derivative of the nucleon correlator in terms of the quark mass 
from the correlation between the nucleon correlator and the scalar density operator (the 
Feynman-Hellmann theorem). Since the staggered fermion has a remnant chiral symmetry, 
there is no problem with the operator mixing. On the other hand, there is a subtlety due 
to the artificial fourth root of the fermion determinant necessary for the staggered fermions, 
for which the Feynman-Hellmann theorem is modified. Their result appears to be slightly 
higher than ours. 
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VII. CONCLUSIONS 



In this paper we calculate the nucleoli strange quark content on the lattice directly from 
the nucleon three-point function in two- flavor QCD. Chiral symmetry is exactly preserved by 
employing the overlap fermion formulation on the lattice. This is crucial in the calculation of 
the strange quark content in order to avoid large contaminations from the operator mixing 
effects, that were missing in many of the previous calculations. 

The lattice calculation of the disconnected diagram is technically challenging. In this 
work we attempted various options of the all-to-all propagator technique and the low-mode 
averaging together with the source and sink smearings. By optimizing those, we could finally 
obtain the nonzero signal at each quark mass; the value extrapolated to the physical quark 
masses is away from zero by 1.5 standard deviation. 

The results for pr and y are in good agreement with our previous estimate using the 



spectrum method 



201 ]. and favor small strange quark content y « 0.05, which is an order 



of magnitude smaller than previous lattice calculations without respecting chiral symmetry, 
which we now believe unreliable. 

For more realistic calculations, we must include the dynamical strange quark in the 
simulation. Such a calculation is already underway using both the spectrum and direct 



methods 
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36j. It is also interesting to extend this study to other baryon observables 



involving disconnected quark loops, such as the strange quark spin fraction of the nucleon. 
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